# =============================================================================
# Statistical Analysis Code for:
# "Chronic adaptive versus conventional DBS response patterns in Parkinson's
#  disease: A pilot randomized crossover trial"
#
# Visualization: Violin plots for clinical outcomes comparison
# Author: Jun Tanimura, MD, MSc.
# Created: 2025-03-08
# Last Updated: 2025-08-07
# =============================================================================

library(tidyverse)

# Read CSV data
df <- read.csv("data.csv", stringsAsFactors = FALSE)
if(!"PatientID" %in% names(df)){
  df$PatientID <- paste0("Patient_", 1:nrow(df))
}

# Identify evaluation items from columns starting with "P_"
p_cols <- names(df)[str_detect(names(df), "^P_")]
eval_items <- sub("^P_", "", p_cols)
cat("Extracted Evaluation Items:\n")
print(eval_items)
# Exclude specific items if needed (e.g., HY, UPDRS_tremor_score, UPDRS_rigidity_score)
eval_items <- eval_items[ ! eval_items %in% c("HY", "UPDRS_tremor_score", "UPDRS_rigidity_score") ]

# Specify the order of evaluation items
eval_items_ordered <- c("Mean_corr_ON_time", "Mean_corr_OFF_time", 
                        "Mean_corr_Dys_sev_time", "Mean_corr_Dys_weak_time", "ADL", "UPDRS_part_1", "UPDRS_part_2", "UPDRS_part_3", "UPDRS_part_4", 
                        "UPDRS_total", "MMSE", "PDQ.39")

# Keep only items that exist in the data
eval_items <- eval_items_ordered[eval_items_ordered %in% eval_items]

cat("Final Evaluation Items (in specified order):\n")
print(eval_items)

# Data transformation for pooled analysis
# Pool all patients' aDBS/cDBS data and compare by treatment mode
pooled_data_long <- map_dfr(eval_items, function(item) {
  col_P <- paste0("P_", item)
  col_C <- paste0("C_", item)
  col_A <- paste0("A_", item)
  
  df %>%
    select(PatientID, Group, all_of(c(col_P, col_C, col_A))) %>%
    rowwise() %>%
    mutate(temp = list(
      tibble(
        Treatment = c("Preoperative", "cDBS", "aDBS"),
        Score = c(get(col_P), get(col_C), get(col_A))
      )
    )) %>%
    unnest(cols = c(temp)) %>%
    mutate(Evaluation = item) %>%
    ungroup() %>%
    select(PatientID, Group, Evaluation, Treatment, Score)
})

# Fix factor order
pooled_data_long <- pooled_data_long %>%
  mutate(Treatment = factor(Treatment, levels = c("Preoperative", "cDBS", "aDBS")),
         Group = factor(Group, levels = c("1", "2")),
         Evaluation = factor(Evaluation, levels = eval_items))

# Add fill_color column for box plot
pooled_data_long <- pooled_data_long %>%
  mutate(fill_color = case_when(
    Treatment == "Preoperative" ~ "darkgray",
    Treatment == "cDBS" ~ "steelblue", 
    Treatment == "aDBS" ~ "darkred",
    TRUE ~ "gray"
  ))

# Draw lines only for patients with all three treatment conditions
valid_patients <- pooled_data_long %>%
  group_by(PatientID, Evaluation) %>%
  summarise(n_treatments = n_distinct(Treatment[!is.na(Score)]), .groups = "drop") %>%
  filter(n_treatments == 3) %>%
  pull(PatientID)

# Create pooled violin plot
p_pooled <- ggplot(pooled_data_long, aes(x = Score, y = Treatment)) +
　geom_violin(aes(fill = fill_color), color = "white", alpha = 0.6, scale = "count", width = 0.8, adjust = 3.0 ) +
  #geom_boxplot(aes(fill = fill_color), color = "white", alpha = 0.7, width = 0.8) +
  geom_path(
    data = pooled_data_long %>% filter(PatientID %in% valid_patients),
    aes(group = PatientID),
    color = "gray", alpha = 0.2, linewidth = 0.8
  ) +
  geom_point(aes(color = fill_color), size = 1.5, alpha = 0.8, 
             position = position_jitter(height = 0.1, width = 0)) +
   stat_summary(fun = median, geom = "point", size = 3, color = "black", shape = 18) +
  facet_wrap(vars(Evaluation), scales = "free_x", ncol = 3, strip.position = "top") +
  scale_fill_identity() +
  scale_color_identity() +
  scale_y_discrete(limits = rev(levels(pooled_data_long$Treatment))) +
  labs(title = "Clinical Outcomes Comparison",
       subtitle = "Distribution and individual patient trajectories across treatment conditions",
       x = "Clinical Score", y = "Treatment Condition") +
  theme_classic() +
  theme(
    strip.background = element_blank(),
    strip.text = element_text(face = "bold", size = 10),
    axis.text = element_text(size = 9),
    axis.title = element_text(size = 11, face = "bold"),
    plot.title = element_text(size = 13, face = "bold"),
    plot.subtitle = element_text(size = 10, color = "gray40"),
    panel.grid = element_blank(),
    legend.position = "none"
  )

print(p_pooled)

# Save the plot
ggsave("violin_plot.pdf", 
       plot = p_pooled, device = "pdf", width = 5, height = 4.5 )

# Optional: Create legend explanation
cat("\nLegend:\n")
cat("Gray: Preoperative baseline\n")
cat("Blue: Conventional DBS (cDBS)\n") 
cat("Red: Adaptive DBS (aDBS)\n")
cat("Gray lines: Individual patient trajectories\n")